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It is well-known that the Liouville equation of statistical mechanics is restricted to systems where 
the total number of particles (N) is fixed. In this paper, we show how the Liouville equation 
can be extended to systems where the number of particles can vary, such as in open systems 
or in systems where particles can be annihilated or created. A general conservation equation 
for an arbitrary dynamical variable is derived from the extended Liouville equation following 
Irving and Kirkwood's 2 technique. From the general conservation equation, the particle number 
conservation equation is obtained that includes general terms for the annihilation or creation of 
particles. It is also shown that the grand canonical ensemble distribution function is a particular 
stationary solution of the extended Liouville equation, as required. In general, the extended 
Liouville equation can be used to study nonequilibrium systems where the total number of 
particles can vary. 
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I. Introduction 

Modern molecular theories of equilibrium thermodynamics and nonequilibrium thermody- 
namics are based almost entirely on the Liouville equation which describes the conservation of an 
ensemble of phase points as they move through phase space according to a specified Hamiltonian 1 . 
Unfortunately, the Liouville equation is based on a fixed number of particles and cannot describe 
systems where the total number of particles varies. For example, the Liouville equation can- 
not describe the statistical mechanical, grand canonical ensemble used in the development of 
thermodynamic properties of open, equilibrium systems. In this letter, we extend the Liouville 
equation to systems where the number of particles varies. To our knowledge this has not been 
done before and has restricted the theoretical development of nonequilibrium thermodynamics. 
These systems include open systems, where the particle number can vary because of open bound- 
aries, or closed or open systems where particles can be annihilated or created either classically 
by chemical reactions or due to quantum mechanical effects. 

II. The Extended Liouville Equation 

The Liouville equation describes the behavior of a collection or ensemble of phase points as 
they move through a multidimensional space, or phase-space. Each phase point represents the 
position (or generalized coordinates g^) and momentum (or conjugate momentum pi) of all N 
molecules or particles in the system. The phase points tend to be concentrated in regions of phase 
space where it is most likely to find the N particles with a certain momentum and position. Thus, 
the Liouville density function, pi, can be interpreted (aside from a normalization constant) as a 
probability density function, i.e., pLdq N dp N is proportional to the probability of finding a phase 
point in a multidimensional region between (q^, p^) and (q^ + dq N , p N + dp N ) at any time, t. 
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Here and throughout we use the short-hand notation q w = q±, ■ • • , q 3 N and p N — pi, - • • ,P3N- It 
can be readily shown that the Liouville density function obeys the conservation equation 1 



0PL(PW,*) | f 



dp L dH dpL dH 



(1) 



_ dqi dpi dpi dqi 

where H (p N , q N , t) is the Hamiltonian or total energy of the iV-particle system. 

It is to be noted that the Liouville equation describes the phase space behavior of the density 
function Pl(p N , <l N , t) under the conditions where the number of particles TV is a constant. To 
extend the Liouville equation to a system where the number of particles or dimensionality can 
vary, consider the extended space (p N ,q N ,N,t) shown in Fig.l. The total members of the 
ensemble in extended space is a constant denoted by M, and the density of phase points is 
denoted by p N {p N , q N , N, t). Thus, we have that 



E / PN(p N ,q N ,N,t)dp N dq N = M (2) 

N 



As shown in Fig.l, the extended space can be divided into "canonical" hyperplanes with the 
number of phase points at time t in each hyperplane denoted by M Nl , M N2 , etc., where 



M = E M Ni (3) 

Ni 

The phase points in each hyperplane have the same number of particles or dimensionality 
and, hence, we call these "canonical" hyperplanes. As the system evolves in time, the phase 
points can move within the plane, if the dimensionality (JV) does not change, or "jump" between 
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planes, if the dimensionality (N) of the phase point changes. In any event, the total number of 
phase points M remains constant with time. 

Now, it is not possible to straightforwardly extend the Liouville equation by considering 
differential changes in q^, and iV (the extended space), since as iV changes so does the 
dimensionality of the phase space coordinates (<l N ,p N ). However, if we define a new function 



p(N, t)=J J p N (p N , q N , N, t)dp N dq N = M N (4) 

then, in order to conserve the total number of phase points M, it can be readily shown (Appendix 
A) that p must obey the one-dimensional conservation equation 



where we have introduced a "velocity" function, v(N,t), 



v(N,t) = v + (N,t) -v~(N,t) (6) 

which physically represents the net fractional number of phase points leaving the N th hyperplane 
per unit time (Appendix A); v + represents the fraction of phase points leaving the N th hyperplane 
to the (N + l) th hyperplane, and v~ are the fraction per unit time leaving to the (N — l) th 
hyperplane. Now, within each hyperplane we introduce the Liouville density p L (p N ,q N ,t; N), 
defined as a conditional density function 
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The phase points that move within the hyperplane over a differential time, dt, obey the usual, 
Liouville equation, Eq.(f), 



dt 




dpL_dB_ _ dp L OH 

dqi dpi dpi dqi 







(8) 



Thus, the space-time behavior of the extended Liouville density function Pn{p n , <l N , N, t) is 
uniquely determined through Eqs.(5), (7) and (8). Physically, we can envision that phase points 
move in each hyperplane over a differential time, dt, according to Eq.(8). At the end of the 
differential time, some phase points "jump" into new canonical hyperplanes according to Eq.(5). 
Thus, the number of phase points in each canonical hyperplane with a fixed dimensionality 
continually change with time. Although both the dimensionality change and time change are 
discrete, the continuous representation given by Eq.(5) should approximately hold for iV large 
and dt — > 0. 

Multiplying Eq.(5) by pi, Eq.(8) by p, and combining using Eq.(7) we can write an extended 
Liouville equation as 



dps (<\ N ,V N ,N,t) 
dt 



+ P(N,t) 



dp L OH dp L OH 

dqi dpi dpi dqi 



+ PL ( q N ,p N ,t;N)^[pv(N,t)} = (9) 



d\np N ^ d\np L dH d\np L dH 

dt Y % ^ d Pi d 1i 



In the next section, we derive a general conservation (transport) equation for an arbitrary 
dynamic variable where the total number of particles can vary. From this we obtain a particle 
number conservation equation which includes a general term accounting for the annihilation or 
creation of particles. 

III. General Conservation Equation for Variable Particle Number Systems 

In Cartesian coordinates, the extended Liouville equation can be written 



dp N (r N ,p N ,N,t) 
dt 



I i 



+p 1 (r\p< v ,t;JV)^[MJV,f)]=0 (11) 



where we will now normalize the density function p by the total ensemble number M (fixed), 



N 



(12) 



J J p L dr N dp N = 1 (13) 



N 



E / / PNdr N dp N = 1 (14) 



Extending Irving and Kirkwood's formalism 2 , the "conservation equations" for variable par- 
ticle number systems can be obtained by first considering an arbitrary quantity a(r N ,p N , N) 
that does not depend explicitly on time. The average or expectation value of a is defined as 



< « E / / «( rJV > p", N)p N dv N dp N (15) 

N 

Now, multiplying the extended Liouville equation by a and integrating over all (r N , p N ) and 
summing over all N yields 



d < a > - f f P i 9a 



dt N 



m J J , m or; 



N 



9a • NjN 



Ep rEPLFi—d^dp 



dpi 



+ E4(^) / fap L dr N dp N (16) 



N dN 



where we have used the conditions 



PL 







VipL 







> as ri,pi 



oo 



(17) 



aViPL 







Equation (16) is the general conservation equation for a. 

Particle Number Conservation 

Again, extending Irving and Kirkwood's formalism 2 for variable particle number systems, 
the equation for particle number conservation (ordinary number density) can be obtained from 
Eq.(16) by setting 



N 




E 





Thus, 




r n 



< a >= n(r, t) = y^p 

N 



N 




N 




(19) 



N 
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where ipi(r,t; N)dr is the conditional probability of finding a particle between r and r + dr, 
given a system of iV total particles at time t, i.e., 



^(r,t; N) = f Pl (r,p,t-N)dp 



(20) 



Note that for a system at equilibrium 



ipi = — , equilibrium 



(21) 



where V is the volume of the system; thus 



I N 
n(r, t) = — ^2 pN = — , equilibrium 



(22) 



where is the equilibrium average number of particles in the volume V. 



Analyzing the remaining terms in Eq.(16) 



N 



Pi da_ 
m <9r, 



N iN 



dr^dp 



d_ 
dr 



N J 171 



d_ 
dr 



[n(r,*)v (r,i)] 



(23) 



where we have defined the number average velocity as 



v (r,t) = 



(24) 



and, for the final term, 



E 

N 



J Jap L dr N dp N 



= E 



N 



N^(r,t; N) 



(25) 



Summarizing, the particle number conservation equation becomes 



M ipv) 



N^(r,t; N)=0 



(26) 



The third term on the left-hand side accounts for the creation or annihilation of particles in the 
system as considered in more detail in the example below. 

Example. First- Order Particle Annihilation 



As a simple example of the application of Eq.(26) consider the problem of particle annihilation 
where the fractional number of particles leaving the N th hyperplane to the (iV — l) th hyperplane 
is a constant k. Thus, Eq.(5) becomes 
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or, in "finite N" form (Appendix A) 



M *- f) -k[p(N + l,t)-p(N,t)] = (28) 



dt 

Consider that initially all phase points are in the N hyperplane, i.e., 



p(iV o ,0) = l (29) 



and 



p(N, 0) = , N e [0, N — 1] (30) 
Around the N hyperplane, the population balance equation is 



+ *«".,*> =0 (31) 



Thus, using the initial state Eq.(29), we obtain 



p(N ,t) = e~ kt (32) 

Now, performing a balance around the N — 1 hyperplane using Eq.(32) leads to 
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dp(N -l,t) 
dt 



+ kp(N - 1,0 = ke 



-kt 



(33) 



and, using the initial condition Eq.(30), we obtain 



p{N - l,t) = kte 



-kt 



(34) 



Continuing this process, for the N hyperplane we must have 



p (iV , t ) = M^ e - fe * 



to! 



(35) 



where m = N — N, m e [0, A?o — 1]. We note that for AT — > oo, 



m=0 



m! 



(36) 



as required. Note that for small or finite values of iVo, we must consider a population balance 
around the N = hyperplane (total particle annihilation plane) in order to satisfy the normal- 
ization condition. 

Now consider the annihilation/creation term in the particle number conservation equation, 
Eq.(26). We have that 



d dp 

dN {pv) = ~m =k 



(N - N) 
kt 



(37) 
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Substituting into the third term on the l.h.s of Eq.(26) leads to the first-order decay kinetic 
expression, under the conditions kt » (N — N), as 



^ + ^'K)+EWi = ° » kt»(N -N) (38) 



N 



or 



^r- + ^-( nv ^ + kn = ' kt»(N -N) (39) 



Note that for small times compared to (N — N)/k the decay kinetics are not first-order, even 
for a constant rate of phase point hopping, k. This example helps illustrate the generality of the 
extended Liouville equation. Next, we address the question as to whether the extended Liouville 
equation leads to an equilibrium behavior consistent with the grand canonical ensemble. 

IV. Equilibrium Behavior. Grand Canonical Ensemble 

Under equilibrium conditions, dp/dt = 0, the general stationary solution to Eq.(5) can be 
written as 



p(N,t) = C 1 f(N) (40) 



v(N,t)=C 2 [f(N)r 1 (41) 

where C\ and C 2 are constants to be determined and f(N) is an arbitrary function of N. Now, 
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it is already well-known that a particular stationary solution to Eq.(8), under the conditions 
dpi/dt = 0, can be written in dimensionally correct form as 3 



Pl(p N , q"; N) = ^exp {-PH(p N , q N )} (42) 



where h is Planck's constant that has the units of action (length x momentum), H(p N ,q N ) is 
the Hamiltonian or total energy, and C 3 and (3 are constants to be determined. 

The constant (5 can be determined by requiring that the average kinetic energy on each 
hyperplane be a uniform constant (canonical ensemble on the hyperplane), i.e., in Cartesian 
coordinates 



(l/2m)/^p L (r",p")dp" _ 3 

" 2 T (43) 



and 



37V -i 9 



H{v\ V N ) = Y,\ P - + ^ N ) (44) 



where $(r ) is the total potential. Thus, we obtain 



fi = w (45) 



Combining Eqs.(40) and (42) gives the extended equilibrium density function as (C 4 = CiC 3 ) 
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p N (p N ,<l N ,N) = ^exp{-/^(p w ,q^)}/(iV) (46) 



The constant C 4 can be determined from Eq.(2) as 



r = (A7) 

4 E N ^f(N)JJe W {-(3H(P N ^ N )}dp N d^ { ] 



Finally, normalizing p^ with the total number of phase points M leads to the well-known form 
for the (classical) probability distribution function for the grand canonical ensemble as 5 , 



p N (p N , q N , N) = [E^Nlh^fWexp {- H ^ T qN) } (48) 



where the (classical) grand partition function E^ rand is given by 



E Grand = Y^ E ™ n f(N) (49) 



N 



where E^ n is the classical canonical partition function 



E% n = [Nlh 31 ^]- 1 J J exp j-/3if (p", q") } dp N dq N (50) 



Now, it is well-known from most probable distribution methods 6 that 



/(TV) = e 
15 



(51) 



for the equilibrium grand canonical ensemble. All that we can state from the extended Liouville 
equation is that Eq.(51) is a particular stationary solution to Eq.(5) 4 . We are presented with the 
same difficulty as in the equilibrium solution to the ordinary Liouville equation, Eq.(29), i.e., the 
lack of proof of the uniqueness of Eq.(29) under equilibrium conditions - a well-known problem 
in statistical mechanics. 

Finally, we note that to complete the grand canonical ensemble description, we must still 
specify A. This procedure is already well-known and is generally accomplished by comparison to 
the thermodynamic relation 6 . 



TdS = dU + pdV - /jdN 



(52) 



where \i is the chemical potential 




(53) 



which leads to 




(54) 



and an expression for the entropy, S, in the grand canonical ensemble as 



S=-~ 



U Nfi 



+ k\nE' 



(55) 



T T 
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where U and N are the average internal energy and particle number over the grand canonical 
ensemble 



The results and procedures given in Eqs.(52) - (57) are already well established, and we 
present them only for completeness. 

IV. Conclusions 

The importance of the analysis given here lies in the evolution equations, Eqs.(5) - (8), or 
Eqs.(9) or (10), that describe, in general, the nonequilibrium behavior of the extended Liouville 
density p^(p N ,q N , N,t) for systems where the number of particles is not fixed. A general con- 
servation equation can be derived, following Irving and Kirkwood's paradigm 2 , for systems where 
the total number of particles can vary. As an example, the general particle number conservation 
equation was obtained that was shown to include a general term for particle annihilation or 
creation. It was also demonstrated that the grand canonical ensemble distribution is a partic- 
ular stationary solution to the extended Liouville equation. As with the equilibrium solution 
to the ordinary Liouville equation, uniqueness of the variable particle number solution under 
equilibrium conditions remains to be demonstrated. It is to be noted that the general station- 
ary solution to the extended Liouville equation, Eqs.(40) and (41), demonstrates that v(N) is 
smallest where p is largest, i.e., near the most probable state. Conversely, phase points move 




A' 



(56) 




N 



(57) 
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or change their dimensionality more rapidly from regions where p is small, i.e., near the least 
probable state. Further treatment of the extended Liouville equation in nonequilibrium systems 
is under current investigation. 
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Appendix A. Conservation Equation for p(N,t) 

Consider a simple population balance around the N th hyperplane (see Fig.l). Since phase 
points cannot be created or destroyed, we must have 



dp(N,t) 
dt 



+ [p(N, t)v + (N, t) - p(N - 1, t)v + (N - 1, t)} 
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-[p(N + 1, t)v~(N + l,t)- p(N, t)v~(N, t)} = (Al) 

where v + (N,t) represents the fraction of phase points in the N th hyperplane leaving to the (iV + 
l) th plane per unit time, v~(N,t) represents the fraction of phase points in the N th hyperplane 
leaving to the (N — l) th plane per unit time, etc. Now, using a finite difference representation of 
a derivative, we can write 



d(pv+) ^ p(N, t)v + (N, t) - p{N - 1, t)v+(N - 1, t) 
dN N-(N-l) 

Thus, Eq.(Al) becomes in the limit of large N 



or, 



where 



(A2) 



dp d(pv + ) d{pv ) _ 



v = v — v (A5) 
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This figure "Figurel.GIF" is available in "GIF" format from: 



http://arXiv.org/ps/physics/9809039v2 



